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1.  ABSTRACT 

A  laser-illuminating  imaging  system  operating  in  the  presence  of  atmospheric  turbulence  will  encounter  several 
sources  of  noise  and  diffraction  induced  errors  that  can  make  postprocessing  the  images  difficult.  As  the  beam 
propagates,  turbulence  induced  tilt  will  cause  the  beam  to  wander  off  axis  from  the  target.  This  is  especially 
troublesome  in  imaging  satellites  due  to  most  turbulence  being  closer  to  the  Earth’s  surface  and  greatly  affecting  the 
beam  in  the  early  stages  of  propagation.  Additionally,  the  returning  beam  convoluted  with  the  target  will  encounter 
turbulence  induced  tilt  that  appears  as  apparent  movement  of  the  target  between  image  frames.  This  results  in 
varying  beam  intensities  at  the  target  surface  between  imaging  frames  that  can  affect  registration  algorithms.  Using 
simulated  data,  an  algorithm  using  expectation  maximization  and  least  squares  techniques  was  developed  that  has 
the  ability  to  estimate  both  the  tilt  of  the  pulsed  laser  beam  and  the  apparent  movement  of  the  object  between 
incoherent  frames  separately  and  produce  a  superior  estimate  of  the  target.  The  results  from  this  algorithm  can  be 
used  to  reduce  the  effects  of  beam  wander  and  thus  increase  the  SNR  of  postprocessed  images. 

2.  INTRODUCTION 

When  high  resolution  imagery  is  desired  from  laser  radar  (LADAR)  imaging  systems  there  are  several  factors  that 
can  limit  the  systems  performance  such  as  diffraction  from  the  LADAR  optics,  atmospheric  turbulence  and  laser 
beam  speckle.  These  factors  can  severely  distort  the  image  quality  and  reduce  the  resolution  of  the  measured  data  in 
each  frame  as  shown  in  Fig.  1.  Due  to  operating  conditions  and  factors  such  as  cost,  size  and  weight,  an  adaptive 
optics  approach  may  not  be  feasible  for  all  situations.  This  is  where  the  benefits  of  a  postprocessing  algorithm  can 
be  exploited  to  improve  the  quality  of  LADAR  imagery.  The  algorithms  studied  in  this  report  can  provide  near  real 
time  estimates  of  the  target  image  with  the  distortions  such  as  speckle,  blurring  and  defocus  significantly  reduced 
when  multiple  frames  of  data  are  available. 

In  a  LADAR  application,  diffraction  due  to  the  atmospheric  turbulence  results  in  primarily  tilt,  blur  and 
defocus  which  account  for  over  93%  of  phase  error  [1].  Atmospheric  turbulence  causes  random  time  delays  in  light 
as  it  propagates  through  the  atmosphere,  using  Fourier  optics  this  time  delay  or  tilt  in  the  propagation  field  can  be 
represented  as  a  spatial  shift  in  the  image  field  [2].  Each  pulse  of  the  laser  beam  is  randomly  shifted  to  a  different 
position  on  the  target  and  the  returning  field  after  propagation  is  again  shifted  and  blurred.  Additionally,  in  a 
LADAR  system,  speckle  is  a  significant  source  of  noise.  Speckle  is  caused  by  the  coherency  of  the  illuminating 
laser  source,  combined  with  the  rough  surface  of  the  target  [3].  Each  frame  of  data  will  contain  independent 
intensity  fluctuations  that  appear  as  bright  and  dark  spots  as  a  result  of  laser  speckle. 
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3.  PREVIOUS  WORK 


Previously,  a  blind  deconvolution  algorithm  has  been  developed  to  improve  the  quality  of  imagery  affected  by  the 
atmosphere’s  optical  transfer  function  (OTF)  and  laser  speckle  by  time  averaging  multiple  deconvolved  and 
registered  frames  [4].  This  algorithm  uses  the  Ayers-Dainty  blind  deconvolution  technique  [5],  however  its 
algorithm  doesn’t  take  into  account  the  beam  shifting  off  axis  between  each  frame  of  data.  Another  technique 
commonly  used  in  blind  deconvolution  problems  is  the  expectation  maximization  (EM)  algorithm  developed  by 
Dempster,  Laird  and  Rubin  [6].  The  EM  algorithm  is  an  iterative  approach  to  computing  the  maximum  likelihood 
estimate  of  the  unknown  variables.  A  benefit  when  working  with  EM  algorithms  is  that  their  convergence  is  assured 
since  the  algorithm  is  guaranteed  to  increase  the  likelihood  function  at  each  iteration. 

This  paper  uses  the  EM  algorithm  approach  to  develop  a  set  of  equations  that  can  iteratively  solve  for  the 
complete  data  by  solving  for  the  target  or  complete  data,  the  global  shifts  and  the  beam  shifts  separately  in  each 
frame. 


4.  SYSTEM  ANALYSIS 

The  imaging  system  simulated  in  this  report  represents  a  LADAR  system  that  uses  a  coherent  laser  source.  The 
target  is  illuminated  by  the  laser  beam  and  the  returning  field  is  distorted  by  atmospheric  turbulence  causing  tilting, 
blurring  and  other  higher  order  [4]  effects  on  the  returned  image.  Additionally  the  intensity  of  the  field  will  contain 
random  speckle  fluctuations.  The  returned  image  is  limited  by  the  size  of  the  beam  and  thus  beam  wander  results  in 
each  frame  being  illuminated  at  a  different  position  and  with  a  different  intensity  than  the  previous  frame. 

An  expression  for  the  image  obtained  from  the  system  is  shown  in  Eq.  1  where  ik  is  the  k,h  measured  data 
frame  which  is  the  true  image,  o,  multiplied  by  the  beam,  b,  and  convoluted  with  h,  the  atmospheric  point  spread 
function  (PSF).  The  global  shifts  are  represented  in  the  PSF  as  ak  and  /?fc  and  the  beam  shifts  are  represented  as  yk 
and  sk. 


ik  =  2_,  2_,  o(z,  w)b(z  —  yk,w  —  sk)h(x  -  z  -  ak,y  -  w  -  ft) 
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4.1  Atmospheric  Turbulence 


Turbulence  in  Earth’s  atmosphere  is  caused  by  random  variations  in  temperature  and  air  motion  which  changes  the 
refractive  index  of  the  air  [7].  As  optical  waves  propagate  in  Earth’s  atmosphere,  the  wave  is  distorted  by  the 
changes  in  the  refractive  index  of  the  air  it  is  traveling  through.  The  variance  in  the  phase  of  the  field  [1]  as  a  result 
of  the  turbulence  is  described  by  Eq.  2  where  D  is  the  diameter  of  the  receiver  aperture  and  r0  is  the  coherence 
diameter.  The  coherence  diameter  or  seeing  parameter  is  typically  used  to  describe  the  optical  quality  of  the 
atmosphere  and  is  typically  around  5-10  cm  for  averaging  viewing  sites  and  up  to  20  cm  for  the  best  viewing  sites. 
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The  phase  variance  due  to  tilt  [1]  in  one  axis  is  given  by  Eq.  3,  this  variance  is  doubled  when  looking  at  both  axis. 
Nearly  87%  of  the  total  phase  variance  is  a  result  of  tilt  with  image  distorting  effects  such  as  blurring  and  defocus 
making  up  the  rest.  Tilt  can  be  successfully  removed  by  an  image  registration  algorithm.  However,  the  higher  order 
image  distortions  with  a  residual  phase  variance  given  in  Eq.  4  are  compensated  for  using  a  blind  deconvolution 
algorithm. 


4.2  Noise 
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The  majority  of  noise  in  a  LADAR  system  can  be  attributed  to  one  of  the  following  sources:  photon  counting  noise, 
laser  speckle  and  background  noise.  Laser  speckle  is  the  result  of  imaging  a  surface,  rough  on  the  scale  of  an  optical 


wavelength,  using  the  coherent  light  from  a  laser  [8],  The  rough  surface  of  an  object  cause  incoming  light  to  reflect 
off  at  random  different  angles  resulting  in  a  multitude  of  amplitude  spread  functions  produced  from  each  point  on 
the  target  interfering  with  one  another  at  the  detector.  The  variance  of  the  laser  speckle,  shown  in  Eq.  5,  has  a 
negative  binomial  distribution  [8]  and  is  dependent  of  the  coherency  of  the  light  and  the  expected  number  of  photons 
received.  An  example  of  this  phenomenon  is  shown  in  Fig.  2b,  the  image  has  a  random  intensity  modulation  that 
severely  distorts  the  image  quality  when  compared  to  the  image  obtained  from  an  incoherent  light  source  shown  in 
Fig.  2a. 

^speckle  =  E[Nsignal}  (l  +  E^nal^j  5 

Where, 

•  E[Nsignal  ]is  the  expected  number  of  photons 

•  Mis  the  coherency  parameter  of  the  light,  1  =  fully  coherent,  D  =  fully  incoherent 


(a)  (b) 

Fig.  1  Incoherent  imaging  and  coherent  imaging  with  only  speckle  noise 


Speckle  noise  can  be  reduced  by  using  a  time  average  of  properly  registered  images.  Depending  on  the  coherence  of 
the  light,  the  speckle  pattern  introduced  to  each  image  can  be  especially  troublesome  when  registering  the  images. 
If  image  quality  is  poor,  the  algorithm  might  not  properly  register  each  image  frame  and  will  blur  the  image. 
Flowever,  a  blind  deconvolution  algorithm  used  to  remove  the  effects  of  the  atmosphere  should  be  able  to  improve 
the  quality  because  the  blurring  effect  from  improperly  registered  frames  is  similar  to  the  blurring  effect  from  the 
atmosphere. 

Photon  counting  noise  occurs  due  to  the  randomness  of  counting  photo-electrons  at  the  detector  and  follows 
a  Poisson  distribution.  Additionally  background  noise  is  the  result  of  radiation  sources  besides  the  illuminating 
beam  that  are  counted  by  the  detector.  The  amount  of  background  noise  can  be  measured  by  taking  images  with  the 
illuminating  laser  not  on. 

5.  METHOD 

The  proposed  approach  to  reduce  the  effects  of  the  atmospheric  turbulence  and  speckle  pattern  of  imagery  obtained 
from  a  LADAR  system  is  to  develop  an  EM  algorithm  that  accounts  for  the  global  shift  and  the  beam  shift  in  each 
frame.  This  paper  describes  the  mathematics  in  developing  this  algorithm  and  then  the  implementation  issues  with 
the  mathematical  solution.  Going  into  the  derivation  it  is  assumed  that  the  beam  shape  is  known  and  can  be 
described  with  a  Gaussian  intensity  profile.  Additionally  the  shape  of  the  PSF  is  known,  however  the  tilt  or  shifts 
are  not  and  also  the  background  radiation  can  be  measured.  These  assumptions  are  not  unrealistic  for  many  LADAR 
applications. 

5.1  Multiframe  Postprocessing  Algorithm 

Similar  to  many  deconvolution  algorithms,  the  expectation  maximization  blind  deconvolution  algorithm  is  modeled 
after  Poisson  statistics  due  to  the  ease  of  working  with  Poisson  statistics.  Following  the  steps  outline  by  Dempster, 
Liard  and  Rubin  [6],  first,  a  statistical  model  for  the  incomplete  data  is  defined  in  Eq.  6.  It  is  known  that  the 
incomplete  data,  d,  is  an  image  array  of  independent  Poisson  random  variables  multiplied  by  the  beam  with 


unknown  shifts  and  convoluted  with  the  shifted  PSF.  The  measured  background  radiation  is  represented  as  B  in  the 
equation. 


E[d(x,  y)]  =  m  o(z,w)b(z  ■ 


Yk’w  ~  £k)h(x  ■ 


ak,y-w-pk)  +  B 
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The  complete  or  desired  data  is  related  to  the  incomplete  data  through  the  relationship  shown  below  in  Eq.  7. 

dk{x,y)=^Y^dk{x,y,z,w)  +  dB  7 


The  complete  data  is  also  considered  to  be  Poisson  so  that  a  statistical  model  of  the  incomplete  data  can  be  related  to 
the  complete  data  through  Eq.  8  and  9.  Choosing  the  complete  data  to  be  Poisson  is  acceptable  because  the  sum  of 
Poisson  random  variables  is  also  a  Poisson  random  variable  [3]. 

E[dk{x,y,z,w )]  =  o(z,w)b(z  —  Yk’w  ~  £k)h{x  —  z  —  ak,y  —  w  —  ft)  8 


E[dg]  =  B  9 

The  algorithm  is  derived  from  the  probability  mass  function  (PMF)  for  a  Poisson  process  [3]  given  in  Eq.  10  where  k 
is  the  number  of  occurrences  and  1  is  the  rate  of  those  occurrences. 


ka-X 


P[k]  = 


XKe 


k\ 
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Using  Eq.  10  and  applying  the  model  for  the  complete  data  in  Eq.  8,  results  in  the  complete  data  likelihood 
expression: 

P[cTk(x,y,z,w )] 

o(z,  w)b(z  -yk,W  -  £k)h(x  -  Z  -  ak,y  -  W  -  pk)dk(x.y^,w)e-o(z,w)b(z-yk,w-ek)h(x-z-ak,y-w-l}k)d^.y^  1 1 

d^ix.y,  z,  w)! 


Expanding  Eq.  1 1  by  solving  for  all  pixel  points  in  the  image  array  results  in  Eq.  12: 


P  [ dk  ( x ,  y,z,w)  V  x,y,z,w] 


-nnnnn 

k  x  y  z  w 


o(z,  w)b(z  -yk,w-  sk)h(x  -z-akJy-w-  pkyk^,y,z,w)e-oiz,wMz-rioW-ekmx-z-ak,y-w-pk^k(’‘^ 

Tk(x,y,z,w)\ 
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Following  the  EM  steps,  the  natural  log  is  taken  to  get  the  log-likelihood  in  Eq.  13.  Greatly  simplifying  this 
expression  is  that  the  products  are  turned  into  summations  when  the  natural  log  is  taken. 


L(.o,ak,pk,Yk,£k )  =  ms  dk{x,y,z,w )  ln(o(z,  w)h(z  —  w  —  £k)h(x  —  z  —  ak,y  —  w  —  /?k)) 

k  x  y  z  w 

—  o(z,w)b(z  —  Yk’w  ~  £k)h{x  —  z  —  ak,y  —  w  —  /?k)  —  In  ( dk{x,y,z ,  w))  ! 
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The  expectation  step  of  the  EM  algorithm  takes  the  conditional  expectation  of  the  complete  data  log  likelihood  when 
given  the  incomplete  data  and  past  estimates  of  the  unknown  parameters.  This  derivation  is  extensive,  however  the 
resulting  expression  is  given  below: 

o,ak,/3k,Yk,£k)\o 

mm- 

k  x  y  z  w 

-  £k)h{x  -  z  -  ak,y  -w  -  ft))  -  o(z,w)b(z  -Yk.w-  £k)h(x  -  z  -  ak,y  -  w  -  ft) 


ld,dk{x,y),a°kld,ttld,Y°kld. 


cold 

£k 


(z,  w)b(z  —  Ykl 


)] 


1)/i(x  —  z  —  akd,y  —  w  —  pg^dtx.y) 


W) 


In  (o(z,w)b(z  —  Ykiw  14 


Where, 


ikold  =  ^  ^  oold(z,w)b(z  -  Ykld’w  ~  £fcid)ft(x  -  z  -  akld,y  —  w  —  /?£id) 

Z  W 


The  maximization  step  of  the  EM  algorithm  maximizes  the  conditional  expected  value  of  the  complete  data  log- 
likelihood  from  Eq.  14  with  respect  to  the  parameters  being  estimated  for  a  single  frame,  k0.  This  process  is 
completed  for  the  three  different  set  of  unknowns:  the  global  shifts,  the  beam  shifts,  and  the  true  object.  In  each 
instance,  the  terms  not  dependent  on  the  parameters  being  estimated  can  be  dropped.  Again,  due  to  extensive 
calculations  required  to  prove  each  likelihood  expression,  the  results  for  each  of  the  three  unknown  parameter  sets 
are  given  below. 


L ■«*  =  Z  Z  °0M(z' w)&(z  “  yC-  w  Z  Z  T^3)  h(x~z~  <■  y~w~  /CXK*  - z  -  “fee  y-w-ft*o) 

z  w  x  y 

Lr ,  =  Z  Z  °OUiz’  W)^Z  “  <“■ w  -  <“) ln  ( b(z  -  Yk°' W  -  £*°))  Z  Z 1  old(r  vl  ■ 

z  w  xy  k«  y,y> 

-  z  Z  o<-z’  w-)fe^z  _  Yk°’ w  ~ 


h(x~z~  <■ y  -  w  -  PC) 


o(z0,w0)  = 


Efc  Ooid(z0,  w0)h(z0  -  Ykld.  Wo  -  £kld )  Sy  doidXr’y\h(x  -  z0  -  akld.  y  -  w0  -  Kld) 
_ ik  yx,y) _ 

Zfch(z0  -yfc,w0  -£k) 
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These  derived  equations  show  that  a  mathematical  solution  to  solving  for  each  of  the  parameters  individually  is 
possible  by  an  iterative  process  that  first  estimates  the  object,  o,  then  estimates  the  global  shift  parameters,  a  and  P, 
then  estimates  the  beam  shifts,  y  and  e.  This  EM  algorithm  is  guaranteed  to  converge  to  the  correct  solution. 
However  a  problem  exists  in  the  implementation  of  the  two  likelihood  expressions  describing  the  shifts.  In  both 
instances,  MATLAB  was  unable  to  evaluate  the  natural  log  of  the  shifting  beams  or  PSFs  with  enough  accuracy  to 
correctly  estimate  the  shift.  The  non-linear  properties  of  the  natural  log  function  caused  the  changes  between  the 
beam  shifts  or  the  PSF  shifts  from  each  frame  to  be  smaller  than  its  precision.  This  derivation  does  however  prove 
that  a  more  accurate  estimate  of  the  object  can  be  found  if  the  beam  shifts  can  be  estimated. 

5.2  Least  Squares  Likelihood 

An  iterative  least  squares  likelihood  approach  was  taken  to  solve  for  the  global  and  beam  shifts.  Each  iteration 
moves  the  estimate  closer  to  the  solution.  This  technique  is  accomplished  one  frame  at  a  time  for  first  the  global 
shifts  and  then  for  the  beam  shifts.  Fig.  2  shows  the  flow  of  this  technique  for  estimating  the  beam  shifts,  estimating 
the  PSF  shifts  is  exactly  the  same  so  it  is  not  shown  here. 


Fig.  2  Iterative  least  squared  error  algorithm 


6.  SIMULATION  &  RESULTS 


A  combination  of  the  object  estimation  using  the  EM  algorithm  in  Eq.  17  and  the  shift  estimation  using  the  least 
squares  algorithm  was  developed  in  MATLAB  using  the  system  characteristics  in  Table  1.  The  algorithm  was 
tested  using  simulated  data.  All  shifts  were  random  uncorrelated  Gaussian  variables,  the  beam  shifts  had  a  standard 
deviation  of  9  pixels  and  the  global  shifts  had  a  standard  deviation  of  3  pixels. 

Each  frame  of  the  simulated  data  was  created  by  shifting  the  unit  amplitude  Gaussian  beam  by  its  random 
shift  and  multiplied  by  the  target  frame  and  then  convoluted  with  a  non-shifted  PSF.  Fourier  optics  properties  were 
used  to  propagate  the  source  to  the  image  plane.  The  image  detected  was  then  shifted  to  simulate  the  tilt  component 
of  the  PSF.  Speckle  was  added  to  each  frame  by  applying  a  negative  binomial  random  variable  to  each  pixel  point 
in  the  imaging  plane.  All  convolutions  were  done  by  transforming  into  the  Fourier  domain  and  multiplying. 


Table  1  Simulation  parameters 


Receiver  aperture  diameter,  D 

.3  m 

Beam  width,  w0 

.08  m 

Wavelength 

1.55  pm 

Propagation  distance 

10  km 

Coherence  factor,  M 

40 

Number  of  frames 

30 

Max  number  of  photons  in  return  pulse 

100  photons 

Background  light 

1  photon 

Number  of  iterations 

20 

Pixel  shift  resolution 

.2  pixels 

The  root  mean  squared  error  (RMSE)  was  calculated  at  each  iteration  since  simulated  data  used  and  the  truth  image 
was  known  using  Eq.  18,  where  o  is  the  estimated  image  produced  at  that  iteration  and  o  is  the  known  true  image. 


RMSE  = 
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The  RMSE  was  first  calculated  using  the  algorithm  developed  and  then  calculated  again  using  the  same  algorithm 
except  the  beam  shifts  were  not  estimated.  The  results,  in  Fig.  3  and  Table  2,  show  that  when  the  beam  position  was 
estimated,  the  algorithm  did  a  better  job  at  estimating  the  object.  The  estimated  image  at  the  iteration  with  the 
minimum  error  is  shown  in  Fig.  4. 
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Fig.  3  Algorithm  RMSE  at  each  iteration 


Table  2  Minimum  photon  error  with  and  without  beam  estimation 


#  of 

iterations 

Min  RMSE 

With  beam  estimate 

11 

12.26  photons 

Without  beam  estimate 

13 

15.01  photons 

Fig.  4  Estimated  object  at  iteration  1 1 


Since  simulated  data  was  used,  the  actual  beam  and  global  shifts  were  compared  to  their  estimated  values.  The 
mean  of  the  minimum  error  difference  is  shown  in  Table  3.  Interestingly,  the  difference  between  the  actual  and 
estimated  shift  values  is  not  minimized  at  the  same  iteration  as  the  object  but  several  iterations  afterwards.  This 
seems  to  imply  that  the  shift  estimators  performance  improves  as  the  object  estimate  improves. 

Table  3  Error  in  shift  estimation 


Mean  Difference 
Error 

a 

.5939  pixels 

P 

.4137  pixels 

Y 

1.2586  pixels 

8 

.7887  pixels 

7.  FUTURE  WORKS 

There  are  several  processes  in  this  algorithm  that  need  to  be  further  studied.  First,  the  EM  algorithm  originally 
developed  would  be  the  ideal  solution  since  it  is  guaranteed  to  converge  and  can’t  fall  into  a  local  minimum  value 
like  the  least  squares  approach.  This  might  produce  shift  estimation  results  that  are  more  accurate  and  likely  further 
increase  the  performance  of  the  algorithm.  Additionally,  the  iterative  approach  needs  defined  stopping  criteria  when 
working  with  non-simulated  data.  This  might  include  iterating  the  beam  and  global  shifts  until  a  minimum  error  is 
estimated  and  then  going  back  and  again  estimating  the  object  with  these  better  shift  parameters.  Lastly,  noise  in 
this  simulation  was  limited  to  speckle  and  background  noise  and  a  more  accurate  model  should  be  used  to  get  a 
more  accurate  performance  baseline. 


8.  CONCLUSION 

When  imaging  with  a  LADAR  system,  image  quality  is  severely  reduced  due  to  atmospheric  turbulence  and  laser 
speckle  noise.  Previous  algorithms  have  been  able  to  compensate  for  these  degradations  using  blind  deconvolution 
and  image  registration  techniques.  However,  these  algorithms  do  not  take  into  account  the  movement  of  the  beam 
between  frames  due  to  atmospheric  induced  tilt.  The  algorithm  developed  in  this  paper  estimated  the  beam  shifts 
along  with  the  global  shifting  between  frames  to  produce  an  estimate  of  the  object  with  significantly  less  effects 
from  atmospheric  turbulence  and  speckle  noise.  Additionally,  it  was  able  to  estimate  the  beam  and  global  shifts 
while  producing  an  estimated  image  that  had  less  error  than  an  estimate  of  the  object  when  the  beam  was  not 
tracked. 
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